{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "21da874c-7ca8-4111-bdae-9eff71fed10b",
   "metadata": {},
   "source": [
    "# Quantum Sine and Cosine Transforms\n",
    "\n",
    "The quantum Sine and Cosine transforms functions are the quantum analog for the discrete Sine and Cosine transforms. The **unitary** versions of the type I and type II transforms are defined as follows: \n",
    "\n",
    "$$\n",
    "{\\rm DCT}^{(1)}_{jk}(N) = \\alpha_{jk}\\sqrt{\\frac{2}{N-1}} \\cos\\left(\\frac{\\pi j k}{N-1}\\right), \\qquad \n",
    "\\alpha_{jk} = \\left\\{ \\begin{array}{l l}\n",
    "\\frac{1}{\\sqrt{2}} & j = 0,N-1 ,\\\\\n",
    "\\frac{1}{\\sqrt{2}} & k = 0,N-1 ,\\\\\n",
    "1 & \\text{else}\n",
    "\\end{array}\n",
    "\\right.,\n",
    "\\qquad j,k = 0\\dots,N-1 \n",
    "$$\n",
    "\n",
    "$$\n",
    "{\\rm DST}^{(1)}_{jk}(N) = \\sqrt{\\frac{2}{N+1}} \\sin\\left(\\frac{\\pi j k}{N+1}\\right), \\qquad j,k = 0\\dots,N-1 \n",
    "$$\n",
    "\n",
    "$$\n",
    "{\\rm DCT}^{(2)}_{jk}(N) = \\alpha_{jk}\\sqrt{\\frac{2}{N}} \\cos\\left(\\frac{\\pi (j+1/2) k }{N}\\right), \\qquad \n",
    "\\alpha_{jk} = \\left\\{ \\begin{array}{l l}\n",
    "\\frac{1}{\\sqrt{2}} & k = 0 ,\\\\\n",
    "1 & \\text{else}\n",
    "\\end{array}\n",
    "\\right.,\n",
    "\\qquad j,k = 0\\dots,N-1 \n",
    "$$\n",
    "\n",
    "$$\n",
    "{\\rm DST}^{(2)}_{jk}(N) = \\alpha_{jk} \\sqrt{\\frac{2}{N}} \\sin\\left(\\frac{\\pi (j+1/2) (k+1)}{N}\\right), \\qquad \n",
    "\\alpha_{jk} = \\left\\{ \\begin{array}{l l}\n",
    "\\frac{1}{\\sqrt{2}} & k = N-1 ,\\\\\n",
    "1 & \\text{else}\n",
    "\\end{array}\n",
    "\\right.,\n",
    "\\qquad j,k = 0\\dots,N-1 \n",
    "$$\n",
    "\n",
    "The open library includes four functions, following the implementation in Ref. [[1](#QCST)]: "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1ec327b-846f-4a7b-8e92-aa13b48306cd",
   "metadata": {},
   "source": [
    "## QCT and QST of type I\n",
    "\n",
    "Function: `qct_qst_type1`\n",
    "\n",
    "Arguments:\n",
    "\n",
    "- `x`: `QArray[QBit]`\n",
    "\n",
    "The `x` quantum argument is the quantum state on which we apply the transforms, according to the following  unitary on $n\\equiv$`x.len` qubits:\n",
    "\n",
    "$$\n",
    "\\left(\n",
    "\\begin{array}{ccc|c}\n",
    "{} &{} &{} \\\\ \n",
    "  {}&{\\rm DCT}^{(1)}(2^{n-1}+1) & {}& 0\\\\\n",
    "  {} &{} &{} \\\\ \n",
    "  \\hline\n",
    "  {} & 0 & {} & i{\\rm DST}^{(1)}(2^{n-1}-1)\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4aaa84d1-36c5-4eb3-a713-593d9cc509b3",
   "metadata": {},
   "source": [
    "### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bc4aae3d-0cd9-4c55-a068-dc71207c3537",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from classiq import *\n",
    "from classiq.execution import (\n",
    "    ClassiqBackendPreferences,\n",
    "    ClassiqSimulatorBackendNames,\n",
    "    ExecutionPreferences,\n",
    ")\n",
    "\n",
    "NUM_QUBITS = 4\n",
    "execution_preferences = ExecutionPreferences(\n",
    "    num_shots=1,\n",
    "    backend_preferences=ClassiqBackendPreferences(\n",
    "        backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR\n",
    "    ),\n",
    ")\n",
    "\n",
    "\n",
    "np.random.seed(123)\n",
    "cos_data = np.random.rand(2 ** (NUM_QUBITS - 1) + 1)\n",
    "cos_data = cos_data / np.linalg.norm(cos_data)\n",
    "sin_data = np.random.rand(2 ** (NUM_QUBITS - 1) - 1)\n",
    "sin_data = sin_data / np.linalg.norm(sin_data)\n",
    "\n",
    "combined_data = np.append(cos_data / np.sqrt(2), sin_data / np.sqrt(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8ddff41a-dd81-4fb7-9146-5faecc486aaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "@qfunc\n",
    "def main(x: Output[QArray[QBit]]):\n",
    "\n",
    "    prepare_amplitudes(combined_data.tolist(), 0.0, x)\n",
    "    qct_qst_type1(x)\n",
    "\n",
    "\n",
    "qmod = create_model(main, execution_preferences=execution_preferences)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4e9bb384-4af0-4ab4-9e2e-2b8ba243ad5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from classiq import synthesize, write_qmod\n",
    "\n",
    "write_qmod(qmod, \"qct_qst_type1\", decimal_precision=15)\n",
    "qprog = synthesize(qmod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "287db2e3-fff2-4e33-954d-877750de7b60",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = execute(qprog).result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "86a9d52b-de56-454c-af47-feff3dafa641",
   "metadata": {},
   "outputs": [],
   "source": [
    "qct_data = np.zeros(2 ** (NUM_QUBITS - 1) + 1).astype(complex)\n",
    "qst_data = np.zeros(2 ** (NUM_QUBITS - 1) - 1).astype(complex)\n",
    "for sample in res[0].value.parsed_state_vector:\n",
    "    if sample.state[\"x\"] < 2 ** (NUM_QUBITS - 1) + 1:\n",
    "        qct_data[int(sample.state[\"x\"])] += sample.amplitude\n",
    "    else:\n",
    "        qst_data[int(sample.state[\"x\"] - 2 ** (NUM_QUBITS - 1) - 1)] += sample.amplitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f901350d-c5c7-44ee-9aaa-a676af01c2c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dct1(n):\n",
    "    dct = np.array(\n",
    "        [\n",
    "            [\n",
    "                np.cos(np.pi * j * k / (n - 1))\n",
    "                * (np.sqrt(1 / 2) if j == 0 or j == n - 1 else 1)\n",
    "                * (np.sqrt(1 / 2) if k == 0 or k == n - 1 else 1)\n",
    "                for j in range(n)\n",
    "            ]\n",
    "            for k in range(n)\n",
    "        ]\n",
    "    ) / np.sqrt((n - 1) / 2)\n",
    "    return dct\n",
    "\n",
    "\n",
    "def dst1(n):\n",
    "    dst = np.array(\n",
    "        [\n",
    "            [np.sin(np.pi * (j + 1) * (k + 1) / (n + 1)) for j in range(n)]\n",
    "            for k in range(n)\n",
    "        ]\n",
    "    ) / np.sqrt((n + 1) / 2)\n",
    "\n",
    "    return dst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f3b5895c-f181-45ed-a2b4-9ef0049c0682",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "measured result: [ 0.64936034 -0.13662084  0.02159456  0.08089956  0.06722088  0.18669631\n",
      "  0.02254746 -0.01198615  0.1123771 ]\n",
      "expected result: [ 0.64936034 -0.13662084  0.02159456  0.08089956  0.06722088  0.18669631\n",
      "  0.02254746 -0.01198615  0.1123771 ]\n"
     ]
    }
   ],
   "source": [
    "global_phase = np.exp(1j * np.angle(qct_data[0]))\n",
    "measured_cos_res = np.real(qct_data / global_phase)\n",
    "expected_cos_res = (dct1(2 ** (NUM_QUBITS - 1) + 1) @ cos_data) / np.sqrt(2)\n",
    "print(\"measured result:\", measured_cos_res)\n",
    "print(\"expected result:\", expected_cos_res)\n",
    "\n",
    "assert np.allclose(measured_cos_res, expected_cos_res, atol=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6889bc0e-9271-4a23-a776-577da8b7d20e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "measured result: [ 0.57556991  0.04712138  0.22433696 -0.27513447  0.17796799  0.07685908\n",
      "  0.05378552]\n",
      "expected result: [ 0.57556991  0.04712138  0.22433696 -0.27513447  0.17796799  0.07685908\n",
      "  0.05378552]\n"
     ]
    }
   ],
   "source": [
    "global_phase = np.exp(1j * np.angle(qst_data[0]))\n",
    "measured_sin_res = np.real(qst_data / global_phase)\n",
    "expected_sin_res = (dst1(2 ** (NUM_QUBITS - 1) - 1) @ sin_data) / np.sqrt(2)\n",
    "print(\"measured result:\", measured_sin_res)\n",
    "print(\"expected result:\", expected_sin_res)\n",
    "\n",
    "assert np.allclose(measured_sin_res, expected_sin_res, atol=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff351ef8-2860-41cc-9f05-622715607b50",
   "metadata": {},
   "source": [
    "## QCT and QST of type II\n",
    "\n",
    "Function: `qct_qst_type2`\n",
    "\n",
    "Arguments:\n",
    "\n",
    "- `x`: `QArray[QBit]`,\n",
    "- `q`: `QBit`\n",
    "\n",
    "The `x` quantum argument is the quantum state on which we apply the transforms, whereas the single `q` qubit indicates the block, according to the following unitary on $n+1\\equiv$ `x.len` $+1$ qubits:\n",
    "\n",
    "$$\n",
    "\\left(\n",
    "\\begin{array}{c|c}\n",
    "  {\\rm DCT}^{(2)}(2^{n-1}) & 0\\\\\n",
    "  \\hline\n",
    "  0 & -{\\rm DST}^{(2)}(2^{n-1})\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "Function: `qct_type2`\n",
    "\n",
    "Arguments:\n",
    "\n",
    "- `x`: `QArray[QBit]`: the quantum state on which we apply ${\\rm DCT}^{(2)}$.\n",
    "\n",
    "\n",
    "Function: `qst_type2`\n",
    "\n",
    "Arguments:\n",
    "\n",
    "- `x`: `QArray[QBit]`: the quantum state on which we apply ${\\rm DST}^{(2)}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c2f605a-208f-476b-9505-5037d18c5b72",
   "metadata": {},
   "source": [
    "### Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "83f82380-de82-4726-8547-3436ebf5f12b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from classiq import (\n",
    "    Output,\n",
    "    QArray,\n",
    "    QBit,\n",
    "    create_model,\n",
    "    execute,\n",
    "    invert,\n",
    "    prepare_amplitudes,\n",
    "    qct_qst_type2,\n",
    "    qfunc,\n",
    ")\n",
    "from classiq.execution import (\n",
    "    ClassiqBackendPreferences,\n",
    "    ClassiqSimulatorBackendNames,\n",
    "    ExecutionPreferences,\n",
    ")\n",
    "\n",
    "NUM_QUBITS = 4\n",
    "cos_sin_data = np.random.rand(2 ** (NUM_QUBITS - 1))\n",
    "cos_sin_data = cos_sin_data / np.linalg.norm(cos_sin_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bb78bf55-f17e-47f2-ba86-045fb1dec7bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "@qfunc\n",
    "def main(x: Output[QArray[QBit]], q: Output[QBit]):\n",
    "\n",
    "    prepare_amplitudes(cos_sin_data.tolist(), 0.0, x)\n",
    "    allocate(1, q)\n",
    "    H(q)\n",
    "    qct_qst_type2(x, q)\n",
    "\n",
    "\n",
    "qmod = create_model(main, execution_preferences=execution_preferences)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "fd3f5d49-87cc-4dfe-a64d-623e7140f387",
   "metadata": {},
   "outputs": [],
   "source": [
    "write_qmod(qmod, \"qct_qst_type2\", decimal_precision=15)\n",
    "qprog = synthesize(qmod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a13e0b0b-755d-403d-9888-36fd2f9e8ec8",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = execute(qprog).result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0de4ad7b-1152-4e0d-b777-077d57a73c81",
   "metadata": {},
   "outputs": [],
   "source": [
    "qct_data = np.zeros(2 ** (NUM_QUBITS - 1)).astype(complex)\n",
    "qst_data = np.zeros(2 ** (NUM_QUBITS - 1)).astype(complex)\n",
    "for sample in res[0].value.parsed_state_vector:\n",
    "    if sample.state[\"q\"] == 0:\n",
    "        qct_data[int(sample.state[\"x\"])] += sample.amplitude\n",
    "    else:\n",
    "        qst_data[int(sample.state[\"x\"])] += sample.amplitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c0ebc3e1-54a1-47a6-b352-b9a6edacf448",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dct2(n):\n",
    "    dct = np.array(\n",
    "        [\n",
    "            [\n",
    "                np.cos(np.pi * j * (k + 1 / 2) / n) * (np.sqrt(1 / 2) if j == 0 else 1)\n",
    "                for j in range(n)\n",
    "            ]\n",
    "            for k in range(n)\n",
    "        ]\n",
    "    ) / np.sqrt(n / 2)\n",
    "    return dct.T\n",
    "\n",
    "\n",
    "def dst2(n):\n",
    "    dst = np.array(\n",
    "        [\n",
    "            [\n",
    "                np.sin(np.pi * (j + 1) * (k + 1 / 2) / n)\n",
    "                * (np.sqrt(1 / 2) if j == n - 1 else 1)\n",
    "                for j in range(n)\n",
    "            ]\n",
    "            for k in range(n)\n",
    "        ]\n",
    "    ) / np.sqrt(n / 2)\n",
    "    return dst.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "59b5a302-3dc0-446f-846d-b7793f530891",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "measured result: [ 0.65104654 -0.23305325 -0.11473439  0.02595719 -0.04930425  0.03323495\n",
      "  0.06553173  0.01252819]\n",
      "expected result: [ 0.65104654 -0.23305325 -0.11473439  0.02595719 -0.04930425  0.03323495\n",
      "  0.06553173  0.01252819]\n"
     ]
    }
   ],
   "source": [
    "global_phase = np.exp(1j * np.angle(qct_data[0]))\n",
    "measured_cos_res = np.real(qct_data / global_phase)\n",
    "expected_cos_res = (dct2(2 ** (NUM_QUBITS - 1)) @ cos_sin_data) / np.sqrt(2)\n",
    "print(\"measured result:\", measured_cos_res)\n",
    "print(\"expected result:\", expected_cos_res)\n",
    "\n",
    "assert np.allclose(measured_cos_res, expected_cos_res, atol=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "9b110435-b340-4515-8d7d-a255bc59cdf6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "measured result: [ 0.63981132 -0.2180174   0.13530848 -0.08552637  0.02796936 -0.03450769\n",
      "  0.12370004 -0.01455965]\n",
      "expected result: [ 0.63981132 -0.2180174   0.13530848 -0.08552637  0.02796936 -0.03450769\n",
      "  0.12370004 -0.01455965]\n"
     ]
    }
   ],
   "source": [
    "global_phase = np.exp(1j * np.angle(qst_data[0]))\n",
    "measured_sin_res = np.real(qst_data / global_phase)\n",
    "expected_sin_res = (dst2(2 ** (NUM_QUBITS - 1)) @ cos_sin_data) / np.sqrt(2)\n",
    "print(\"measured result:\", measured_sin_res)\n",
    "print(\"expected result:\", expected_sin_res)\n",
    "\n",
    "assert np.allclose(measured_sin_res, expected_sin_res, atol=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2f0a410d-5bf2-4ab7-bd69-3ec7dac1305c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "measured result: [ 0.92071884 -0.32958707 -0.16225893  0.03670901 -0.06972674  0.04700132\n",
      "  0.09267587  0.01771754]\n",
      "expected result: [ 0.92071884 -0.32958707 -0.16225893  0.03670901 -0.06972674  0.04700132\n",
      "  0.09267587  0.01771754]\n"
     ]
    }
   ],
   "source": [
    "@qfunc\n",
    "def main(x: Output[QArray[QBit]]):\n",
    "\n",
    "    prepare_amplitudes(cos_sin_data.tolist(), 0.0, x)\n",
    "    qct_type2(x)\n",
    "\n",
    "\n",
    "qmod = create_model(main, execution_preferences=execution_preferences)\n",
    "\n",
    "write_qmod(qmod, \"qct_type2\", decimal_precision=15)\n",
    "qprog = synthesize(qmod)\n",
    "res = execute(qprog).result()\n",
    "qct_data = np.zeros(2 ** (NUM_QUBITS - 1)).astype(complex)\n",
    "\n",
    "for sample in res[0].value.parsed_state_vector:\n",
    "    qct_data[int(sample.state[\"x\"])] += sample.amplitude\n",
    "\n",
    "global_phase = np.exp(1j * np.angle(qct_data[0]))\n",
    "measured_cos_res = np.real(qct_data / global_phase)\n",
    "expected_cos_res = dct2(2 ** (NUM_QUBITS - 1)) @ cos_sin_data\n",
    "print(\"measured result:\", measured_cos_res)\n",
    "print(\"expected result:\", expected_cos_res)\n",
    "\n",
    "assert np.allclose(measured_cos_res, expected_cos_res, atol=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "ae387709-e7b5-472c-806c-0b712a298b68",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "measured result: [ 0.90482984 -0.30832317  0.19135509 -0.12095255  0.03955464 -0.04880124\n",
      "  0.17493827 -0.02059046]\n",
      "expected result: [ 0.90482984 -0.30832317  0.19135509 -0.12095255  0.03955464 -0.04880124\n",
      "  0.17493827 -0.02059046]\n"
     ]
    }
   ],
   "source": [
    "@qfunc\n",
    "def main(x: Output[QArray[QBit]]):\n",
    "\n",
    "    prepare_amplitudes(cos_sin_data.tolist(), 0.0, x)\n",
    "    qst_type2(x)\n",
    "\n",
    "\n",
    "qmod = create_model(main, execution_preferences=execution_preferences)\n",
    "\n",
    "write_qmod(qmod, \"qst_type2\", decimal_precision=15)\n",
    "qprog = synthesize(qmod)\n",
    "res = execute(qprog).result()\n",
    "qst_data = np.zeros(2 ** (NUM_QUBITS - 1)).astype(complex)\n",
    "\n",
    "for sample in res[0].value.parsed_state_vector:\n",
    "    qst_data[int(sample.state[\"x\"])] += sample.amplitude\n",
    "\n",
    "global_phase = np.exp(1j * np.angle(qst_data[0]))\n",
    "measured_sin_res = np.real(qst_data / global_phase)\n",
    "expected_sin_res = dst2(2 ** (NUM_QUBITS - 1)) @ cos_sin_data\n",
    "print(\"measured result:\", measured_sin_res)\n",
    "print(\"expected result:\", expected_sin_res)\n",
    "\n",
    "assert np.allclose(measured_sin_res, expected_sin_res, atol=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8aeffc6e-526d-4680-a4a8-34c5b534db2b",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "<a id='QCST'>[1]</a>: [Klappenecker, A., & Rotteler M., \"Discrete Cosine Transforms on Quantum Computers\".](https://arxiv.org/abs/quant-ph/0111038)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
